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Chemical thermalization in relativistic heavy ion collisions 
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We compute by numerical integration of the Dirac equation the number of quark-antiquark pairs 
initially produced in the classical color fields of colliding ultrarelativistic nuclei. While the number 
of qq pairs is parametrically suppressed in the coupling constant, we find that in this classical field 
model their production rate is comparable to the thermal ratio of gluons/pairs = 9-/Vf/32. After 
isotropization one thus would have a quark-gluon plasma in chemical equilibrium. 
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Much attention has recently been given to the ther- 
malization properties of QCD matter formed in ultra- 
relativistic heavy ion collisions. Since 1983 |l( hydrody- 
namic analysis have assumed rapid initial thermalization 
and essentially entropy conserving expansion thereafter. 
Now RHIC experiments 0, strongly suggest that this 
is the correct picture. Theoretically, it is straightforward 
to understand the formation of an initial gluonic system, 
but the problem has been its isotropization in momen- 
tum space 0. Weak coupling methods fail if they only 
include collective effects by screening [f| but show great 
promise if they include collective effects caused by the 
anisotropy of the system 0, @ ■ However, this still 
leaves open the chemical equilibration of the system [ljj . 



The physics involved here is as follows. Two large nu- 
clei of radius Ra moving along the light cone in opposite 
directions (Fig. ^) collide with zero impact parameter. 
They correspond to an ensemble of color currents moving 
along the light cone in opposite directions. The stochas- 
tic properties of the transverse color density p a (xp) have 
been extensively studied 0, 0] . The magnitude of the 
charging is described by a parameter combination g/i. 
For given element of the charge ensemble, color fields 
can be computed by solving the Yang-Mills equations 
[D I j i ,F p ' u ] = with initial conditions given by continu- 
ity across the light cone 01 ■ F° r a dilute system (small 
,',,;? . ^ one can compute the gluon fields analytically 



In general, the color fields can be obtained by 
a numerical computation [lil flfij | and then interpreted 
in terms of gluon production. The color fields are space 
and time dependent and a definite quantum mechani- 
cal probability of producing qq pairs is associated with 
them. These are pairs which are nonperturbatively pro- 
duced during the first instants of the collision and which 
will contribute to chemical equilibration. We will com- 
pute their number and various distributions for numeri- 
cally computed gluonic fields. Again, for weak fields this 
computation can be carried out analytically |l7j . 

Assume now that the SU(3) color fields 




FIG. 1: Domains of different time dependences. The fermion 
amplitude is a sum of two terms: one with interaction first 
with the left moving nucleus, then evolution over a region in 
which the gluon field is a gauge transform of the vacuum 
and finally evolution in a nontrivial gluon field Ai,A v . The 
other term is symmetric. 



Ai(r, xy), Ajjir, xt) 1 are known for < t < l/Q s ~ 
0.2 fm/c (RHIC value). A computation of pair produc- 
tion requires a careful discussion of in and out vacua 
and their relation. For the present field configuration 
this has been carried out in The computation, 

which is set up so that one obtains the average number 
of pairs, which is just the quantity we want, proceeds as 



1 Our longitudinal variables are t, z, x^ = {t± z) / 'y/2 = re ±,? Iv2, 
r = Vt 2 - z 2 = V2x+x~, n=\ la(x + /x~), ds 2 = At 2 - dz 2 = 
2dx+dx~ = dr 2 - r 2 dn 2 so that A T = A T = (tA° - zA a )/t = 
(x+A~ + x~A+)/t and A v = -t 2 A^ = zA° - tA 3 = x+A~ - 
x~ A + . We work in the Schwinger gauge A T = 0. 
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follows. Choose an antiquark of momentum q and mass 
m and solve as a function of time the Dirac equation 
with this color field for the spinor ?/> q (i,x) which in 
the distant past is given by the negative energy spinor 
ijjq(t — > — oo, x) = e tq ' x v(q). The time integration brings 
in positive energy components and consists of three 
qualitatively different domains, see Fig. ^ The region 
< is trivial. The regions marked A^\A^ can be 
dealt with analytically |l9( and one obtains an initial 
condition for ?/»q( T = 0, z,x.t) along the positive light 
cones. This rather complicated initial condition, given 
explicitly in Eq. (16) of 0, depends on the Wilson 
lines {/(i)(xt), J7(2)(xt) corresponding to the gauge 
fields of the nuclei, the initial color field Aj(0,xr) and 
on y qi (\T,z,x.T- The spinor ?/> q (i, x) at r > is then 
computed by solving the Dirac equation in the given 
color field forward in time. Finally, one chooses a quark 
momentum p and forms the overlap between a positive 
energy spinor 2 <f) p (x) = e~ lp ' x u(p) and the outcome of 
the time evolution of the negative energy spinor in the 
distant past: 



M T (p,q) = 



rdzd 2 xy 
\Jt 2 + z 2 



4(T,x) 7 V^q(T,x). (1) 



The overlap is computed at fixed r, hence the use of 7 1 ", 
7°7 T = cosh 77 — 7°7 3 sinh?7 = exp(— ?77 7 3 ). This is also 
the reason for the Jacobian factor rdz / \Jt 2 + z 2 in the 
longitudinal integration. We evaluate Eq. (JTJ in the 2- 
dimensional Coulomb gauge diAi = 0. This is the gauge 
condition used in the Abelian case HI and also the one 
used to evaluate the number of gluons in the background 
field. Eq. gives us 



dN 
dy 



dy p d 2 p T dy g d 2 q T 
2 (2tt) 3 2 (2tt) 3 



S(y-y p )\M T (p,q)\\ (2) 



the number of quarks of one flavor of mass m per unit 
rapidity (since an equal number of antiquarks are pro- 
duced, we refer to this quantity as the "number of pairs" 
below). Since the gluon fields are 77-independent, dN/dy 
is independent of y. We shall compute J2J f° r au r but it 
is only after the "formation time" r > l/y^J + rn 2 that 
the produced antiquarks can reinteract. Since one ex- 
pects qr ~ <7 2 M, this limit for light quarks is r > l/(g 2 [i). 

The parameters of the computation are the coupling g 
(constant in this semiclassical set-up; we use the phe- 
nomenologically relevant value g = 2, a s = 0.3) the 
source density parameter [i (depends on atomic number 
A and collision energy \J~i) the nuclear radius Ra and the 
quark mass m (like with g there is nothing in this semi- 
classical set-up which would make m scale dependent). 



2 Whether it is justified to use a free spinor at a finite r in the 
presence of the external field merits further study. 
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FIG. 2: Dependence on proper time r of the number of pairs 
per unit rapidity dN/dy for g 2 fi — 2 GeV and for values of 
quark mass marked on the figure. The lowest curve corre- 
sponds to g 2 fi = 1 GeV. 
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FIG. 3: Dependence of the number of quark pairs on quark 
mass at a fixed proper time, r = 0.25 fm, and for two values 
of pV 



Two relevant parameter combinations are <? 2 /i and 
g 2 fiR,A- The first one, <? 2 /i, is the dominant transverse 
momentum scale of the classical background field, ft 
is related to the saturation scale Q s ; for < Q s the 
gluonic system becomes so dense that nonlinear interac- 
tions limit the growth of its density, numerically in one 
phenomenological model [H Q s « 0.2 GeV^ 0128 y/f' 19 , 
Q s w 1 GeV at RHIC energies and w 2 GeV at LHC en- 
ergies. The dimensionless diluteness parameter g 2 [iRa 
determines the importance of nonlinear strong field ef- 
fects. 

The numerical computation is done on a N%,N L lattice 
so that the total transverse area is (Nto) 2 = 7r(6.7 fm) 2 , 
i.e., the transverse lattice spacing is a = 11.8fm/iVy = 
60/Nt ■ 1/GeV. The results presented in this letter have 
been obtained with N T = 180, N L = 400, dz = 0.2a and 
dr = 0.02a. At each site one has for each color a spinor 
with 4 complex components, i.e., (N c = 3)x2x4x4 = 96 
bytes in single precision, giving a total of 96 ■ 180 2 • 400 = 
1.2 GB. This illustrates the memory requirement of the 
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FIG. 4: Dependence of the number of quark pairs on g 2 n at 
a fixed proper time, r = 0.25 fm, and for quark mass m — 
0.3 GeV. 



calculation. The number of timesteps in the integration 
in order to reach r = 0.25 fm is of the order of 500. 
The numerical method can be tested by varying these 
parameters. Another check of the numerical method is 
to study how for a zero external field the result of zero 
quark pairs is obtained in a rather nontrivial way by the 
contributions of the two paths in FigQ] cancelling each 
other. 

Note that on the transverse lattice one has to use lat- 
tice momenta 



A 2 

=1 



2 



a- 



sin , k = y sin aki, (3) 

2 or f— ' 



for bosons and fermions, aki = 2imi/Nx with —Nt/2 + 
1 < n,- < Nt/2. For fermions the modes iVr/4 < |n^| < 
Nt/2 are doubler modes, which we leave out both from 
the initial condition (q modes) and the projection to pos- 
itive energy states {p modes). For fermions one can ef- 
fectively only go up to k = %/2/a « %/2iV T /60 GeV « 
4 GeV on a 180 2 -lattice. 

Results for the number of pairs are shown in Figs. [21 El 
First, Fig. |2 shows how the pair production amplitude 
M r depends on proper time r for g 2 fi — 2 GeV and for 
a number of quark mass values, for g 2 ^ = 1 GeV the r 
dependence for m = 0.3 GeV is given. A striking feature 
of the result is the instantaneous formation of the pairs 
at small r, followed by gradual increase. To put this in 
perspective, note that in the QED case the production 
amplitude Q can be computed analytically for all r > 
and its square is constant and equal to that of the pair 
production amplitude in |l8j . In Abelian electrodynam- 
ics the production thus takes place instantaneously at 
t = 0, Here there is also a slow increase thereafter and 
one will enter in a domain in which the backreaction of 
qq on the gluon fields should be included. In comparison, 
for g 2 fi — 2 GeV the gluon number grows rapidly to the 
value 0.29ng 2 n 2 R 2 A « 1000 by the time l/g 2 ^ « 0.1 fm 
[T5L flif and saturates thereafter. 
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FIG. 5: Transverse momentum spectrum of (anti)quarks for 
g 2 fi = 2 GeV at a fixed proper time, r = 0.25 fm, and for 
different quark masses. 
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FIG. 6: Transverse momentum spectrum of (anti)quarks for 
quark mass m = 0.3 GeV and for different g 2 fi at a fixed 
proper time, r = 0.25 fm. 



The dependence of the result on quark mass m and on 
g 2 \i is shown in Figs. |3and01 both at r = 0.25 fm. The 
dependence on the quark mass m is surprisingly weak 
compared to the perturbative result. For numerical rea- 
sons we can not go down to quark masses m < 100 MeV. 
One should also note that with the small transverse lat- 
tice that we have used here one cannot study very large 
masses. The mass dependence is stronger for smaller g 2 /j,, 
which is expected. 

The computation of the gluon fields 0, 0] is based 
on the Hamiltonian formalism and thus gives directly the 
transverse energy Ex- Obtaining the multiplicity is based 
on assuming free field dispersion relation. In the case of 
quarks, the pair multiplicity comes directly, and one must 
explicitly compute the transverse momentum spectra in 
order to obtain E^ lls . These are obtained from |J2J by 
fixing y q — and qr and integrating over y p , pt- Results 
for various quark masses at fixed g 2 [i = 2 GeV are shown 
in Fig. |31and for various g 2 ^, at fixed m = 0.3 GeV in 

Fig. El 

One expects gr-spectra to become flatter and smaller 
with increasing quark mass m and flatter and larger with 
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increasing g 2 fx. These qualitative features are seen in the 
result, but especially the mass dependence is rather weak. 

It has conventionally been assumed that the initial 
state of a heavy ion collision is dominated by gluons. This 
is the result e.g. when both quarks and gluons are pro- 
duced in 2 —>■ 2 collisions of collinear partons [2(J • In the 
color glass condensate picture the initial state is entirely 
gluonic and pair production is suppressed by a power of 
g. Our result suggests that quarks could be present in 
comparable numbers. This is understandable since the 
numerical value of the coupling constant g is, in fact, not 
small. Note also that the kinematics are different: the 
calculation of gluon production in this approach reduces 
in the perturbative limit to a 2 — > 1 process, whereas the 
weak field limit of our present computation is a 2 — > 2 
process. To be able to compare quark and gluon produc- 
tion to the same order in g one must compute the first 
correction to the gluon production result, which has not 
yet been done in the color glass condensate picture. A 
first step in this direction is relaxing the assumption of 
boost invariance in solving the gauge field equations of 
motion |2l| . 

The qualitative phenomenological implications of our 
result are as follows. Experiments at RHIC have ob- 
served about 600 charged (~ 1000 charged and neutral) 
particles per unit rapidity. The assumption of entropy 
conservation, which is supported by the success of ideal 
hydrodynamics calculations, then implies as many parti- 
cles per unit rapidity also in the initial state, which we 
are discussing. One has normally assumed that all the 
initial particles are gluons. In this framework, g 2 ji = 2 
GeV would lead to 1000 initial gluons. Now we see from 
Fig. |21 that, for 3 flavors, about 600 pairs would be pro- 
duced. This implies that the initial assumption of gluon 
dominance should be questioned. In chemical thermal 
equilibrium the quark/gluon ratio is 9N{/32 w 1 mean- 
ing that 1000 particles in a unit of rapidity should consist 
of 330 gluons and 330 pairs. These numbers are obtained 
for g 2 /i as 1.2, i.e., the model gives a consistent fit with a 
reduced value of <? 2 /i, smaller saturation scale. 

There is now experimental evidence suggesting that 
QCD matter formed in ultrarelativistic heavy ion colli- 
sions is from the very beginning in local kinetic thermal 
equilibrium. This is what one always has assumed in ap- 
plications of hydrodynamics, but its theoretical proof and 
understanding has been lacking as the process clearly is 
nonperturbative. We have in this letter shown that the 
classical gluon field model with an ensemble of initial con- 
ditions also produces an abundance of quark-antiquark 
pairs. In fact, this number is close to the one dictated by 
chemical equilibrium between quarks and gluons. This 
suggests that the QCD matter formed initially in heavy 
ion collisions could be in full chemical and kinetic thermal 
equilibrium. Experimentally, this would make it more 
probable to observe thermal dilepton radiation in the 
2 GeV mass range at LHC energies. It would be very 



interesting to find also other experimental tests for ini- 
tial chemical equilibration. 
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